ARL-TR-8270  •  JAN  2018 


ARL 

US  Army  Research  Laboratory 


An  Automated  Energy  Detection  Algorithm 
Based  on  Morphological  Filter  Processing  with 
a  Modified  Watershed  Transform 


by  Kwok  F  Tom 


Approved  for  public  release;  distribution  is  unlimited. 


NOTICES 

Disclaimers 

The  findings  in  this  report  are  not  to  be  construed  as  an  official  Department  of  the 
Army  position  unless  so  designated  by  other  authorized  documents. 

Citation  of  manufacturer’s  or  trade  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  thereof. 

Destroy  this  report  when  it  is  no  longer  needed.  Do  not  return  it  to  the  originator. 


ARL-TR-8270  •  Jan  2017 


ARL 

US  Army  Research  Laboratory 


An  Automated  Energy  Detection  Algorithm 
Based  on  Morphological  Filter  Processing  with 
a  Modified  Watershed  Transform 


by  Kwok  F  Tom 

Sensors  and  Electron  Devices  Directorate ;  ARL 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  the  collection  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  the 
burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

This  report  is  the  application  of  morphological  image  techniques  to  the  energy  detection  scenario  of  signals  in  the  RF 
spectrum  domain.  The  algorithm  automatically  establishes  a  detection  threshold  based  on  the  data  of  the  RF  spectrum 
measurement.  Initially,  a  modified  watershed  transform  was  applied  to  form  an  envelope  of  the  spectral  data.  Then  a 
morphological  “opening”  operation  was  applied  to  the  spectral  envelope  data  to  establish  the  noise  floor.  Using  this  noise 
floor  estimation,  an  energy  detection  threshold  was  determined. 


15.  SUBJECT  TERMS 

RF  spectrum,  morphological  image  processing,  detection  threshold  algorithm,  opening,  watershed  transform 

17.  LIMITATION  18.  NUMBER 
OF  OF 

ABSTRACT  PAGES 

UU  74 

Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 

ii 


19a.  NAME  OF  RESPONSIBLE  PERSON 

Kwok  F  Tom 

19b.  TELEPHONE  NUMBER  (Include  area  code) 

(301)  394-2612 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

Unclassified 

Unclassified 

Unclassified 

1.  REPORT  DATE  (DD-MM-YYYY) 

2.  REPORT  TYPE 

January  2018 

Technical  Report 

4.  TITLE  AND  SUBTITLE 

An  Automated  Energy  Detection  Algorithm  Based  on  Morphological  Filter 
Processing  with  a  Modified  Watershed  Transform 


6.  AUTHOR(S) 

Kwok  F  Tom 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

US  Army  Research  Laboratory 

Sensors  and  Electron  Devices  Directorate  (ATTN:  RDRL-SER-E) 
2800  Powder  Mill  Road 
Adelphi,  MD  20783-1138 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED  (From  -  To) 

1  October  2016-30  September  2017 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION  REPORT  NUMBER 

ARL-TR-8270 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT  NUMBER(S) 


Contents 


List  of  Figures  v 

List  of  Tables  v 

Preface  vi 

1.  Introduction  1 

2.  Data  Collection  and  Statistical  Summary  1 

3.  Statistical  Processing  2 

3.1  Statistical  Analysis  2 

3.1.1  Moments  2 

3.1.2  Mean  3 

3.1.3  Variance  3 

3.1.4  Standard  Deviation  3 

3.1.5  Kurtosis  3 

3.1.6  Maximum  4 

3.1.7  Minimum  4 

3.1.8  Median  4 

3.1.9  Rank  Order  Filter  4 

3.1.10  Crest  Factor  (CF)  5 

3.1.11  Running  Median  5 

3.2  Statistical  Summary  6 

4.  Morphological  Image  Processing  7 

4.1  Dilation,  Erosion,  and  Opening  8 

4.2  Watershed  Transformation  9 

5.  Algorithm  10 

6.  Conclusion  11 


Approved  for  public  release;  distribution  is  unlimited. 

iii 


7.  References 


13 


Appendix  A.  MATLAB  Code  15 

Appendix  B.  Graphs  of  Morphological  Processed  RF  Spectrum  Files  29 

Appendix  C.  Graphs  of  RF  Spectrum  Files  Calculated  Detection 

Threshold  47 

List  of  Symbols,  Abbreviations,  and  Acronyms  64 

Distribution  List  65 


Approved  for  public  release;  distribution  is  unlimited. 


IV 


List  of  Figures 


Fig.  1  Spectral  data  file  with  modified  watershed  transform . 9 

Fig.  2  Zoomed- in  examination  of  spectral  data  file  with  modified  watershed 

transform . 10 

List  of  Tables 

Table  1  Summary  of  the  RF  spectrum  collection . 2 

Table  2  Statistical  analysis  of  the  RF  spectrum  measurements . 7 


Approved  for  public  release;  distribution  is  unlimited. 


V 


Preface 


Energy  detection  in  the  RF  spectrum  is  the  most  basic  technique  for  signal 
detection.  Typically,  this  requires  establishing  an  energy  detection  threshold  based 
on  a  noise-only  condition  (i.e.,  no  signal  present  in  the  RF  spectrum).  An 
investigation  of  morphological  image  processing  techniques  was  successful  in 
automatically  generating  an  energy  detection  threshold  based  on  the  RF  spectrum 
data  without  the  requirement  of  having  RF  spectrum  of  noise  only  condition. 

The  technique  of  “opening”  was  very  successful  in  determining  an  energy  detection 
threshold  for  the  RF  spectrum  with  signal  and  noise-only  environments.  A  second 
energy  threshold  detection  algorithm  was  examined  with  the  addition  of  a  semi¬ 
disk  structure  applied  to  the  RF  spectrum.  During  the  development  of  these 
automatic  energy  detection  threshold  generation  algorithms,  the  watershed 
transform  proved  to  be  a  good  candidate  to  apply  to  the  RF  spectrum  prior  to  the 
application  of  the  opening  technique.  This  is  the  third  of  5  reports  that  detail  the 
energy  detection  techniques  examined  with  the  recorded  RF  spectrum 
measurements. 1-4 


1  Tom  K.  An  automated  energy  detection  algorithm  based  on  morphological  filter  processing  with  a  semi¬ 
disk  structure.  Adelphi  (MD):  Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8271. 

2  Tom  K.  An  automated  energy  detection  algorithm  based  on  morphological  and  statistical  processing 
techniques.  Adelphi  (MD):  Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8272. 

1  Tom  K.  An  automated  energy  detection  algorithm  based  on  kurtosis-histogram  excision.  Adelphi  (MD): 
Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8269. 

4  Tom  K.  An  automated  energy  detection  algorithm  based  on  consecutive  mean  excision.  Adelphi  (MD): 
Army  Research  Laboratory  (US);  2018  Jan.  Report  No.:  ARL-TR-8268. 
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1.  Introduction 


Energy  detection  is  the  simplest  method  of  detecting  a  signal  in  the  frequency 
spectrum.  When  doing  so,  a  comparison  is  made  between  the  frequency  spectral 
component  energy  and  a  detection  threshold  level.  Knowledge  of  the  frequency 
spectrum  is  usually  necessary  to  establish  this  detection  threshold  level.  Various 
methods  can  be  used  to  determine  the  detection  threshold  level,  such  as  establishing 
the  system  noise  statistics  offline  and  setting  the  threshold  for  a  given  probability 
of  detection  versus  probability  of  false  alarm. 

The  goal  of  this  study  was  to  develop  an  algorithm  that  can  analyze  the  frequency 
spectrum  and  establish  a  threshold  value  based  on  the  spectral  components.  The 
automated  processing  employs  techniques  from  image  processing  and  statistical 
analysis.  A  combination  of  techniques  from  these  2  areas  resulted  in  an  algorithm 
for  determining  a  detection  threshold  level.  In  an  earlier  evaluation,  the 
morphological  filter  processing  successfully  determined  the  detection  threshold 
level.  In  this  evaluation,  a  modified  watershed  transformation  is  applied  to  the  RF 
spectral  data  prior  to  the  morphological  filtering  process.  A  reduction  in  false 
alarms  was  obtained. 

2.  Data  Collection  and  Statistical  Summary 

The  local  RF  spectrum  was  measured  in  2013  on  the  rooftop  of  building  204  at  the 
US  Army  Research  Faboratory’s  (ARF’s)  Adelphi  location.  An  Agilent  N9342CN 
spectrum  analyzer  and  a  Discone  Antenna  was  used  to  collect  RF  spectrum  data. 
This  spectrum  analyzer  was  operated  under  the  control  of  a  Fabview  software 
program  to  acquire  and  store  data  with  different  resolution  bandwidth  (RBW)  from 
1  KHz  to  1  MHz. 

These  data  files  represent  various  sizes  of  RF  spectral  coverage  from  10  MHz  up 
to  4  GHz.  The  number  of  data  files  for  each  spectral  band  varied.  The  larger  RBW 
files  were  acquired  over  seconds  of  data  acquisition  time  versus  the  small  RBW 
files.  Small  RBW  provides  fine  spectral  resolution,  but  it  has  impact  on  data 
acquisition  time  for  spectral  coverage  and  data  size.  Depending  on  the  RBW  and 
spectral  coverage,  a  data  file  could  require  a  few  hours  of  acquisition. 

Table  1  summarizes  the  data  collection  measurements  files.  The  RF  spectral  bands 
covered  the  spectrum  from  10  MHz  up  to  4  GHz.  In  general,  various  spectral  bands 
were  measured  with  4  RBW  configuration.  Data  file  size  is  inversely  proportional 
to  the  RBW.  Data  file  size  is  proportional  to  the  spectral  band  coverage.  The  data 
size  varied  from  approximately  1  KSample  to  4  MSample  data  points  per  file. 
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Table  1  Summary  of  the  RF  spectrum  collection 


Spectral  band  coverage 

Number  of  RBW 
measurement 

RBW 

1.1-1. 6  GHz 

4 

1  kHz,  10  kHz,  100  kHz, 

2-3  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

3^1  GHz 

3 

1  MHz 

10  kHz,  100  kHz,  1  MHz 

4-6  GHz 

3 

10  kHz,  100  kHz,  1  MHz 

10  MHz-1  GHz 

4 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-2  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-3  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-4  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

100  MHz-1  GHz 

1 

1  MHz 

100  kHz 

3.  Statistical  Processing 

3.1  Statistical  Analysis 

Statistical  analysis  is  the  mathematical  science  dealing  with  the  analysis  or 
interpretation  of  data.  The  data  analyst  uses  a  few  straightforward  statistical 
techniques  as  a  means  of  summarizing  the  collected  data.  These  statistical 
techniques  are  under  the  area  of  descriptive  statistics,  which  is  a  methodology  to 
condense  the  data  in  quantitative  terms. 

In  commercial  prognostics  and  diagnostic  vibrational  monitoring  applications, 
statistical  techniques  that  are  mainly  used  for  alarm  purposes  in  industrial  plants 
are  the  statistical  moments  of  order  2,  3,  and  4.  The  probability  density  function 
(PDF)  of  the  vibrational  time  series  of  a  good  bearing  has  a  Gaussian  distribution 
(also  known  as  a  normal  distribution),  whereas  a  damaged  bearing  results  in  a  non- 
Gaussian  distribution  with  dominant  tails  because  of  a  relative  increase  in  the 
number  of  high  levels  of  acceleration.  These  techniques  can  be  applied  to  the  RF 
spectral  data  with  a  different  interpretation  of  the  results.1 

3.1.1  Moments 

If  these  moments  are  calculated  about  the  mean,  they  are  called  central  statistical 
moments.  The  first  and  second  moments  are  well  known,  being  the  mean  and  the 
variance,  respectively.  These  are  analogous  to  the  first  and  second  area  moments 
of  inertia  with  the  area  shape  defined  by  the  PDF.  The  third  moment  is  termed 
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skewness  and  the  fourth  moment  is  termed  kurtosis.  The  general  equation  for  the 
order  of  moment  is  as  follows: 


N 


Mt 


=  "  xy- 


1=1 


where  p  is  the  order  of  the  moment, 
N  is  the  number  of  data  value, 
i  is  the  index  of  the  data  value,  and 
x  is  the  mean  value  of  the  data  set. 


3.1.2  Mean 

Mean  is  the  most  common  measure  of  a  statistical  distribution.  In  this  case,  mean 
is  the  arithmetic  average  for  a  set  of  measurements. 


i= 1 


3.1.3  Variance 

Variance  is  a  measure  of  the  dispersion  of  a  waveform  about  its  mean — also  called 
the  second  moment  of  the  measurements. 

N 

ff2 = ■  *)2- 
i= 1 


3.1.4  Standard  Deviation 

Standard  deviation  is  a  measure  of  the  variation  of  a  set  of  data  values.  The  standard 
deviation  is  defined  as  the  square  root  of  the  variance  moment. 


a 


\ 

1 

i= 1 


(x  i  -  x y. 


3.1.5  Kurtosis 

Kurtosis  is  the  fourth  statistical  moment,  normalized  by  the  standard  deviation  to 
the  fourth  power.  It  is  a  measure  of  whether  the  data  are  peaked  or  flat  relative  to  a 
normal  distribution.  The  noise  in  the  RF  spectrum  is  typically  considered  to  have  a 
normal  distribution.  The  normal  distribution  has  a  value  of  3. 
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K 


M4 
a*' 

lull (*■  -  XT 

K  =  — - - 

<74 

N 

k = -^Ysxi  - *y- 

i= 1 


3.1.6  Maximum 

“Max”  is  the  largest  value  of  a  set  of  numbers. 

y  —  max[x(n)]. 


3.1.7  Minimum 

“Min”  is  the  smallest  value  of  a  set  of  numbers. 

y  -  min[x(n)]. 


3.1.8  Median 

The  statistical  median  is  an  order  statistic  that  gives  the  “middle”  value  of  a  set  of 
samples.  Median  is  the  middle  value  of  a  set  of  data  values  that  divides  the  set  into 
2  groups.  Half  the  groups  exist  below  and  half  exist  above  this  value. 

3.1.9  Rank  Order  Filter 

Rank  order  filter  is  a  sorting  process  by  which  a  set  of  numbers  is  ordered  from  the 
smallest  to  the  largest  value.  Rank  order  filtering  is  a  nonlinear  filtering  technique 
that  orders  the  contents  of  a  filter  kernel  and  selects  the  sample  indexed  by  rank 
from  the  magnitude  ordered  samples. 

Rank  order  filtering  can  be  summarized  as  follows: 

y  =  xr=1a^(0  . 

where  x (q,  l  —  1,  •••  ,N  is  the  result  of  sorting  the  data  in  ascending  order.  With  this 
definition,  the  max,  min,  and  median  can  be  obtained  from  a  rank  order  filter  as 
follows2: 

Min: 

y  =  T,i=i  , 

at  —  l  for  i  =  1, 
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0  otherwise. 


Max: 


Median: 


y  =  £f=ia^(0 , 

N 

at  =  l  fori  = 

0  otherwise. 

y  =  Sf=i  a{X(j)  , 

N 

di  =  l  fori  = 

0  otherwise. 


3.1.10  Crest  Factor  (CF) 

Crest  factor  (CF)  is  a  measure  of  a  waveform  showing  the  ratio  of  peak  values  to 
the  effective  value.  In  other  words,  CF  indicates  how  extreme  the  peaks  are  in  a 
waveform. 

Crest  Factor  =  -  - peak _ 

Xrms 

Noise  sources  are  characterized  by  their  CF,  which  is  the  peak  to  average  ratio  of 
the  noise.  In  a  technical  bulletin,  XiTRON  reported  CF  values  between  5  and  7  for 
random  noise.3  For  example,  a  5: 1  CF  of  the  noise  voltage  is  201og5  =  14  dB.  This 
is  a  measure  of  the  quality  of  the  noise  distributions  and  one  way  to  measure  its 
Gaussian  nature.4  For  the  purpose  of  algorithm  development,  the  CF  equation  was 
modified  as  follows: 

Max 

Crest  Factor  —  - - - . 

Median 


3.1.11  Running  Median 

To  define  the  running  median  filter,  let  {x}  be  a  discrete  time  series.  The  running 
median  passes  a  window  over  the  sequence  { jc  }  that  selects,  at  each  instant  m,  an 
odd  number  of  samples  to  comprise  the  observation  vector  x(n).  The  observation 
window  is  typically  symmetric  and  centered  at  n,  resulting  in 

x(n )  =  [x(n  —  Nf),  ■■■  ,x(n),  ,x(n  +  Ni)]T. 


Approved  for  public  release;  distribution  is  unlimited. 

5 


where  N1  may  range  in  value  over  the  nonnegative  integers  and  N  —  2N{  +  1  is 
the  (odd  valued)  window  size.  While  processing  such  noncausal  observation 
vectors  has  traditionally  been  referred  to  as  smoothing,  we  loosen  the  terminology 
somewhat  and  refer  to  the  processing  of  both  causal  and  noncausal  observations  as 
simply  filtering.  The  median  filter  operating  on  the  input  sequence  {x}  produces 
the  output  sequence{y},  where  at  time  index  n 

y(n)  =  MED  [x(n)] 

=  Median  value  of  [x(n  —  Nf),  •••  ,x(n),  •••  ,x(n  +  Nf)]. 

That  is,  the  samples  in  the  observation  window  are  sorted,  and  the  middle,  or 
median,  value  is  taken  as  the  output.5 

3.2  Statistical  Summary 

There  were  31  different  groupings  of  the  RF  spectrum  data  measurements.  The 
number  of  data  files  under  each  of  the  main  groupings  was  varied.  For  the  purpose 
of  developing  the  algorithm,  only  a  single  data  file  was  selected  from  each  group. 
Each  data  file  was  processed  to  obtain  the  following  characteristics:  RBW,  mean, 
standard  deviation,  median,  max,  min,  kurtosis,  and  CF.  A  summary  of  the  results 
is  in  Table  2. 

The  following  results  were  noted: 

•  The  smaller  the  RBW,  the  lower  the  noise  floor.  The  equation  for  thermal 
noise  power  is  P  —  kTB,  where  B  is  bandwidth.  In  this  case,  the  RBW  of  1 
kHz  has  the  lowest  noise  value.  The  RBW  of  1  MHz  has  the  highest  noise 
value. 

.  Each  10-fold  increase  in  bandwidth  results  in  a  10-dB  increase  in  noise 
power.  This  relationship  is  illustrated  in  this  data  set. 

•  The  calculated  mean  and  median  values  are  very  close  for  a  given  RF 
measurement  configuration. 

•  The  CF  for  noise-only  data  files  was  on  the  order  of  approximately  10  to 
13.  Noise-only  data  files  were  estimated  by  visually  inspecting  the  spectrum 
plot. 
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Table  2  Statistical  analysis  of  the  RF  spectrum  measurements 


Filename 

Spectral  Band 

RBW 

Mean 

SD 

Median 

Max 

Min 

Kurtosis 

CF 

Air_test_l.  lGHz_1.6GHzb_03_28_14_06_33_22 

1.1-  1.6  GHz 

1000 

-112.694 

3.8444 

-112.3 

-101.2 

-138.4 

3.5571 

11.1 

Ai  r_test_l.  lGHz_l.  6GHza_03_27_14_07_14_05 

1.1-  1.6  GHz 

10000 

-102.12 

3.7226 

-101.8 

-91.69 

-124.5 

3.5134 

10.11 

Ai  r_test_l.  1G  Hz_l.  6G  Hzc_03_31_14_06_47_26 

1.1-  1.6  GHz 

100000 

-90.9109 

2.597 

-90.75 

-83.13 

-101.2 

3.0217 

7.62 

Ai  r_test_l.  lGHz_1.6GHzd_04_01_14_06_54_03 

1.1-  1.6  GHz 

1000000 

-80.5634 

3.5741 

-80.745 

-55.36 

-89.8 

14.6861 

25.385 

A  i  r_t  e  st_2G  H  z_3G  H  za_06_05_14_07_15_58 

2- 3 GHz 

1000 

-111.219 

4.4991 

-111.1 

-77.55 

-139 

6.8918 

33.55 

A  i  r_t  e  st_2G  H  z_3G  H  z  b_05_29_14_06_28_28 

2- 3 GHz 

10000 

-101.006 

4.1808 

-100.8 

-74.67 

-122 

6.1995 

26.13 

A  i  r_t  e  st_2G  H  z_3G  H  zc_05_29_14_04_09_27 

2- 3 GHz 

100000 

-89.5666 

3.4154 

-89.68 

-64.27 

-100.3 

13.2282 

25.41 

A  i  r_t  e  st_2G  H  z_3G  H  zd_05_28_14_00_06_18 

2- 3 GHz 

1000000 

-79.6014 

3.0753 

-79.6 

-58.74 

-88.88 

8.4466 

20.86 

A  i  r_t  e  st_3G  H  z_4G  H  z  b_06_12_14_06_30_17 

3-  4  GHz 

10000 

-100.51 

3.7072 

-100.2 

-89.68 

-123.7 

3.5901 

10.52 

A  i  r_te  st_3G  Hz_4G  H  zC_06_ll_14_06_44_40 

3  -  4  GHz 

100000 

-89.1625 

2.5513 

-89.075 

-80.68 

-102.3 

3.0774 

8.395 

A  i  r_te  st_3G  Hz_4G  H  zd_06_09_14_07_09_5 1 

3-  4  GHz 

1000000 

-79.0828 

2.445 

-78.94 

-72.76 

-88.5 

2.9651 

6.18 

A  i  r_te  st_4G  Hz_6G  Hza_04_24_14_07_05_26 

4-6GHz 

1000 

-107.559 

4.1793 

-107.2 

-94.31 

-136.3 

3.3696 

12.89 

A  i  r_t  e  st_4G  H  z_6G  H  z  b_04_28_14_06_23_12 

4-  6 GHz 

10000 

-97.245 

4.0398 

-96.95 

-84.05 

-119.7 

3.3249 

12.9 

Ai  r_te  st_4G  Hz_6G  H  zc_04_29_14_06_48_42 

4-6GHz 

100000 

-85.9125 

3.082 

-85.77 

-75.84 

-98.78 

2.9068 

9.93 

A  i  r_te  st_10M  Hz_lG  Hza_04_17_14_07_09_30 

10  MHz  - 1  GHz 

1000 

-110.348 

8.6126 

-112.1 

-28.37 

-140.3 

6.4788 

83.73 

Ai  r_te  st_10M  Hz_lG  Hz  b_04_21_14_07_17_44 

10  MHz  - 1  GHz 

10000 

-100.168 

8.4775 

-102 

-27.95 

-126.7 

6.7602 

74.05 

A  i  r_te  st_10M  Hz_lG  Hzc_04_22_14_06_39_58 

10  MHz  - 1  GHz 

100000 

-87.948 

8.8113 

-90.74 

-27.28 

-103.6 

6.513 

63.46 

A  i  r_te  st_10M  Hz_lG  Hzd_04_23_14_06_36_16 

10  MHz  - 1  GHz 

1000000 

-78.1054 

8.3327 

-80.59 

-27.82 

-89.88 

7.5282 

52.77 

A  i  r_t  e  st_10M  H  z_2G  H  zd_03_06_14_08_52_04 

10  MHz -2  GHz 

1000 

-111.275 

6.9798 

-112.1 

-27.34 

-138.3 

9.5085 

84.76 

A  i  r_te  st_10M  Hz_2G  HzC_02_27_14_06_43_55 

10  MHz -2  GHz 

10000 

-100.894 

7.0157 

-101.8 

-26.58 

-125.9 

10.4224 

75.22 

A  i  r_t  e  st_10M  H  z_2G  H  z  B_02_25_14_06_35_53 

10  MHz -2  GHz 

100000 

-89.242 

6.7487 

-90.59 

-26.09 

-102.4 

12.7961 

64.5 

A  i  r_te  st_10M  Hz_2G  Hz_02_20_14_06_55_36 

10  MHz -2  GHz 

1000000 

-78.8448 

7.3203 

-80.6 

-26.54 

-92.41 

11.3819 

54.06 

A  i  r_t  e  st_10M  H  z_3G  Hz_03_l  1_14_06_32_5 1 

10  MHz -3  GHz 

1000 

-111.373 

6.1356 

-111.8 

-29.42 

-140.5 

11.0474 

82.38 

A  i  r_t  e  st_10M  H  z_3G  H  z  B_03_13_14_06_32_45 

10  MHz -3  GHz 

10000 

-100.787 

6.3591 

-101.4 

-25.36 

-125.9 

11.7254 

76.04 

A  i  r_te  st_10M  Hz_3G  HzC_03_18_14_06_32_08 

10  MHz -3  GHz 

100000 

-89.1614 

6.0257 

-90.14 

-24.16 

-102.4 

15.3376 

65.98 

A  i  r_t  e  st_10M  H  z_3G  H  zd_03_20_14_08_09_07 

10  Mhz  -  3  GHz 

1000000 

-78.6492 

6.7418 

-80.8 

-25.81 

-90.38 

13.0673 

54.99 

Ai  r_te  st_10M  Hz_4G  Hza_04_10_14_07_01_22 

10  MHz -4  GHz 

1000 

-111.158 

5.6813 

-111.4 

-28.22 

-143.1 

10.9573 

83.18 

A  i  r_te  st_10M  Hz_4G  Hz  b_04_ll_14_06_10_18 

10  MHz -4  GHz 

10000 

-100.838 

5.5851 

-101.1 

-26.85 

-128 

12.2576 

74.25 

A  i  r_te  st_10M  Hz_4G  Hzc_04_15_14_07_01_09 

10  MHz -4  GHz 

100000 

-89.4431 

5.1076 

-90.03 

-25.05 

-103.5 

18.317 

64.98 

Ai  r_te  st_10M  Hz_4G  Hzd_04_16_14_06_42_27 

10  MHz -4  GHz 

1000000 

-79.0396 

5.5441 

-79.89 

-26.72 

-91.64 

16.9073 

53.17 

A  i  r_t  e  st_100M  H  z_lG  H  z_02_19_14_07_20_46 

100  MHz  - 1  GHz 

100000 

-88.0772 

9.0701 

-90.82 

-23.91 

-102.2 

6.994 

66.91 

4.  Morphological  Image  Processing 

A  description  of  a  technique  for  automatically  estimating  the  noise  floor  spectrum 
was  given  at  a  conference  back  in  1997. 6  This  technique  is  based  on  applying  the 
morphological  binary  image  processing  operators  to  the  RF  spectrum.  The 
technique  works  well  in  both  flat  and  nonflat  noise  floor  spectra.  Morphology 
image  processing  is  a  set  of  nonlinear  operations.  Ready  et  al.  note  that  “humans 
are  good  at  estimating  the  noise  floor  spectrum  by  ‘eyeballing’  a  spectral  plot. 
Intuitively,  we  separate  the  spectral  humps  from  the  noise  floor  spectrum  by 
eliminating  those  parts  of  the  spectrum  shape  that  are  due  to  signals  and  visually 
draw  in  the  noise  floor  spectrum.”6 

Morphology  is  a  broad  set  of  processing  techniques  that  process  images  based  on 
shapes.  Morphological  operations  apply  a  structuring  element  to  an  input  image, 
creating  an  output  image  of  the  same  size.  In  a  morphological  operation,  the  value 
of  each  pixel  in  the  output  image  is  based  on  a  comparison  of  the  corresponding 
pixel  in  the  input  image  with  its  neighbors.  By  choosing  the  size  and  shape  of  the 
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neighborhood,  you  can  construct  a  morphological  operation  that  is  sensitive  to 
specific  shapes  in  the  input  image. 


4.1  Dilation,  Erosion,  and  Opening 


The  most  basic  morphological  operations  are  dilation  and  erosion.  Dilation  adds 
pixels  to  the  boundaries  of  the  objects  in  an  image,  while  erosion  removes  pixels 
on  object  boundaries.  The  number  of  pixels  added  or  removed  from  the  objects  in 
an  image  depends  on  the  size  and  shape  of  the  structuring  element  used  to  process 
the  image.  In  the  morphological  dilation  and  erosion  operations,  the  state  of  any 
given  pixel  in  the  output  image  is  determined  by  applying  a  rule  to  the 
corresponding  pixel  and  its  neighbors  in  the  input  image.  These  rules  are  known  as 
the  erosion  and  dilation: 


Erosion: 


Dilation: 


AeB=„EBA-‘- 


a@b  =  u  ABa . 

a  G  A  a 


The  opening  of  A  by  B  is  obtained  by  the  dilation  of  A  by  B,  followed  by  erosion  of 
the  resulting  structure  by  B. 


Opening: 


A  o  B  =  (A  ©  S)©B. 


The  application  of  the  opening  technique  on  the  RF  spectrum  is  the  basis  for  noise 
floor  estimation.  The  estimation  of  the  noise  floor  is  used  as  a  reference  level  to  set 
a  threshold  detection  level.7 


In  signal,  statistical,  and  image  processing,  the  minimum  and  maximum  operators 
are  typically  encountered.  Gil  and  Kimmel  note,  “In  mathematical  morphology,  the 
result  of  such  an  operator  is  referred  to  as  the  erosion  (or  dilation)  of  the  signal  with 
a  structuring  element  given  by  a  pulse  of  width  p.”8 

For  the  1-D  case,  this  reduces  to  a  simple  filter  of  just  providing  the  max  or  min 
value  of  a  set  of  values.8 


1-D  max  filter:  Given  a  sequence  x0,  ••• ,  xn_1  and  an  integer  p  >  1,  compute 


Yi  —  max  xi+i 

o <j<p  J 


for  i  =  0 ,  •••  ,n  —  p. 
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1-D  min  filter:  Given  a  sequence  xQ,  ■■■ , xri_1  and  an  integer  p  >  1,  compute 

Vt  =  mm  xi+j, 

o  <J<p 

for  i  =  0,  •••  ,n  —  p. 

4.2  Watershed  Transformation 

In  morphology  image  segmentation,  there  is  a  technique  called  watershed 
transformation.  It  has  various  definitions  and  implementations.  In  concept,  it 
partitions  an  image  into  disjoint  regions  so  that  these  regions  have  the  same 
properties.  Imagine  a  landscape  with  various  hills  and  valleys.  As  the  area  is 
flooded  by  water,  the  low-lying  areas  will  fill  up,  thus  creating  a  uniform  region. 
Roerdink  and  Meijster  note,  “When  the  water  level  has  reached  the  highest  peak  in 
the  landscape,  the  process  is  stopped.  As  a  result,  the  landscape  is  partitioned  into 
regions  or  basins  separated  by  dams,  called  watershed  lines.”9 

There  are  various  implementations  of  the  watershed  transformation.  The  purpose 
of  applying  the  watershed  transformation  was  to  obtain  the  envelope  of  the  spectral 
data  array.  For  the  spectral  data  array,  the  transformation  was  modified  to  a  highly 
localized  implementation.  A  search  is  conducted  for  2  adjacent  peak  amplitudes  in 
the  spectral  data  array.  The  spectral  data  values  are  replaced  between  the  2  adjacent 
peak  locations  with  the  smaller  of  the  2  peak  amplitude  values  if  the  spectral  data 
are  less  than  that  value.  Figure  1  shows  a  spectral  data  file  with  the  modified 
watershed  transform.  The  red  curve  is  the  modified  watershed  transform  curve.  The 
blue  curve  is  the  measured  spectral  data.  Figure  2  shows  a  zoomed-in  portion  of  the 
spectral  data  file. 


rtowl  OOOOOO-Avg-O— r#f  »•  v*t »-30  OOOOOO-J tin  wj lu*-20 


Fig.  1  Spectral  data  file  with  modified  watershed  transform 
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rt>w*1000000-avg«©— f*<_J*v*1«-30.000000--*ltn_valij#«20 


Fig.  2  Zoomed-in  examination  of  spectral  data  file  with  modified  watershed  transform 

5.  Algorithm 

The  algorithm  process  was  developed  in  MATLAB.  Appendix  A  is  the  code  used 
to  process  and  generate  the  enclosed  RF  spectrum  signatures  with  the 
corresponding  results  of  the  morphological  processing  and  resultant  threshold 
detection  level.  The  following  is  a  description  of  the  detection  thresholding  level 
generation  based  on  morphological  filter  processing: 

1)  Determine  the  RBW  of  the  spectral  data. 

2)  Determine  some  statistics  on  the  spectral  data  file:  median,  max,  and  CF. 

3)  Starting  from  the  beginning  of  the  spectral  data  array,  determine  the  first  2 
peak  amplitude  locations  and  values.  Replace  the  data  array  between  these 
2  peak  locations  with  the  smallest  peak  amplitude  value  if  the  data  value  is 
less  than  the  smallest  peak  value.  Search  for  the  next  peak  amplitude 
location  and  value.  Using  this  location  with  the  last  peak  location,  replace 
any  spectral  data  value  in  the  valley  between  the  latest  2  peak  locations  with 
the  smaller  amplitude  of  the  2  peaks.  Continue  to  search  and  replace  until 
the  end  of  the  spectral  data  array. 

4)  Perform  a  sliding  running  median  on  the  envelope  spectral  data  array  for  a 
window  size  of  5.  The  start  and  end  of  the  spectral  data  are  somewhat  an 
issue.  In  this  case,  the  first  5  points  used  the  median  value  determined  at  the 
first  valid  value.  The  last  5  points  used  the  median  value  determined  at  the 
last  valid  point. 
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5)  Perform  an  erosion  operation  on  the  median  filtered  spectral  data  as 
computed  in  step  4.  Starting  from  the  lowest-frequency  component,  create 
the  eroded  spectral  array  with  an  initial  sliding  window  of  size  2.  Basically, 
calculate  the  minimum  value  for  the  given  window  size,  fill  in  that  position, 
and  shift  the  processing  window  one  position  in  frequency  position  to 
calculate  the  next  value.  Continue  on  across  the  entire  median  filtered 
spectral  data  array.  For  the  last  value,  use  the  previous  valid  value  that  is 
calculated. 

6)  Perform  a  dilation  operation  on  the  eroded  spectral  array  as  calculated  in 
step  5.  In  this  case,  the  calculations  are  performed  starting  at  the  highest 
frequency  and  continue  to  the  lowest  frequency.  Dilation  is  simply  the 
maximum  value  in  the  windowed  data  set.  The  initial  window  size  is  2.  In 
this  case,  the  maximum  value  is  determined  from  the  sliding  window  and 
creating  the  dilated  spectral  array. 

7)  Calculate  statistics  on  the  dilated  spectral  array  from  step  6:  mean,  median, 
max,  and  CF. 

8)  Repeat  the  morphological  filtering  of  erosion  and  dilation  on  the  spectral 
array  by  increasing  the  window  size  by  1  for  every  pass.  For  example,  on 
the  second  pass,  the  window  size  is  3  for  the  eroding  and  dilation  operations. 
The  erosion  and  dilation  operations  are  applied  to  the  dilated  spectral  array 
from  step  6.  Repeat  this  process  as  long  as  the  CF  is  greater  than  10  and  the 
mean  power  greater  than  0.01  on  each  pass  of  the  morphological  filtering. 

9)  If  the  CF  is  not  greater  than  10,  then  perform  the  morphological  filtering 
operations  until  the  average  power  is  greater  than  0.01. 

10)  Once  the  convergence  has  been  met,  use  the  resulting  morphological 
filtered  data  array  (i.e.,  the  dilated  spectral  array)  to  form  the  threshold  as 
follows: 

Threshold  —  Morphological  Filtered  Spectral  Array  +  25  — 

„  _  „  _  .  ( RBW\ 

2.9xl01oglo(— ) 

6.  Conclusion 


In  general,  this  enhanced  algorithm  works  well  in  determining  the  detection 
threshold  level. 

Appendix  B  displays  the  graphs  of  the  morphological  processing  for  each  of  the 
selected  RF  spectrum  data  files.  The  red  curve  is  the  result  of  the  morphological 
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processing  overlaid  on  top  of  the  RF  spectrum  signature.  The  processing  results  in 
a  curve  that  converges  to  the  bottom  of  the  noise  components. 

Appendix  C  displays  the  results  of  the  threshold  overlaid  on  top  of  the  RF  spectrum 
signature.  The  graphs  are  intended  to  provide  the  reader  with  a  qualitative  sense  of 
the  effectiveness  of  the  automatic  detection  threshold  generation  algorithm.  False 
alarms  are  reduced  with  the  prefiltering  of  the  modified  watershed  technique 
compared  with  the  morphological  filtering- only  case.  The  actual  processing  time 
and  number  of  iterations  through  the  morphological  filtering  routine  are  highly 
dependent  on  the  RBW  and  spectral  coverage  band.  A  small  RBW  over  a  large 
coverage  band  requires  more  time  to  complete  the  process.  The  threshold  detection 
level  obtained  works  well  for  signals  present  or  noise-only  RF  spectrum  data  files. 
The  red  curve  is  the  threshold  when  added  to  the  morphological  processed  array  as 
shown  in  Appendix  B. 

Overall,  a  visual  inspection  of  the  graphs  shows  that  the  algorithm  works  well.  The 
false  alarm  for  the  10-kHz  RBW  is  higher  than  desired.  The  threshold  parameters 
can  be  modified  to  obtain  the  desired  level  of  detection  versus  the  false  alarm. 
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Appendix  A.  MATLAB  Code 
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function  WatershedMorphologicalFilterDetection4report () 


%  Make  selection  for  data  file  to  process 
%  SELECT  =  1  to  31 

SELECT  =  31; 


dir{l}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  1.1GHz  1.6GHza\'; 

filename{l}  =  'Air  test  1.1GHz  1 . 6GHza  03  27  14  07  14  05'; 

dir{2}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  1 . 1GHz  1 . 6GHzb\ ' ; 

filename{2}  =  'Air  test  1.1GHz  1 . 6GHzb  03  28  14  06  33  22'; 

dir{3}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  1 . 1GHz  1 . 6GHzc\ ' ; 

filename{3}  =  'Air  test  1.1GHz  1.6GHzc  03  31  14  06  47  26'; 

dir{4}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  1 . 1GHz  1 . 6GHzd\ ' ; 

filename{4}  =  'Air  test  1.1GHz  1.6GHzd  04  01  14  06  54  03'; 

dir{5}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  2GHz  3GHza\ ' ; 

filename{5}  =  'Air  test  2GHz  3GHza  06  05  14  07  15  58'; 

dir{6}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  2GHz  3GHzb\ ' ; 

filename{6}  =  'Air  test  2GHz  3GHzb  05  29  14  06  28  28'; 

dir{7}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzc\ ' ; 

filename{7}  =  'Air  test  2GHz  3GHzc  05  29  14  04  09  27'; 

dir{8}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  2GHz  3GHzd\ ' ; 

filename{8}  =  'Air  test  2GHz  3GHzd  05  28  14  00  06  18'; 

dir{9}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air  test  3GHz  4GHzb\ ' ; 

filename{9}  =  'Air  test  3GHz  4GHzb  06  12  14  06  30  17'; 
dir {10}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  3GHz  4GHzC\ ' ; 
filename{10}  =  'Air  test  3GHz  4GHzC  06  11  14  06  44  40'; 
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dir{ll}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  3GHz  4GHzd\ ' ; 
filename{ll}  =  'Air  test  3GHz  4GHzd  06  09  14  07  09  51'; 

dir{12}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  4GHz  6GHza\ ' ; 
filename{12}  =  'Air  test  4GHz  6GHza  04  24  14  07  05  26'; 

dir{13}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  4GHz  6GHzb\ ' ; 
filename{13}  =  'Air  test  4GHz  6GHzb  04  28  14  06  23  12'; 

dir{14}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  4GHz  6GHzc\ ' ; 
filename{14}  =  'Air  test  4GHz  6GHzc  04  29  14  06  48  42'; 

dir{15}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz_lGHza\ ' 
filename{15}  =  'Air  test  10MHz  lGHza  04  17  14  07  09  30'; 

dir{16}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz^lGHzb\ ' 
filename{16}  =  'Air  test  10MHz  lGHzb  04  21  14  07  17  44'; 

dir{17}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz_lGHzc\ ' 
filename{17}  =  'Air  test  10MHz  lGHzc  04  22  14  06  39  58'; 

dir{18}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz_lGHzd\ ' 
filename{18}  =  'Air  test  10MHz  lGHzd  04  23  14  06  36  16'; 

dir{19}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  2GHz\'; 
filename{19}  =  'Air  test  10MHz  2GHz  02  20  14  06  55  36'; 

dir{20}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz_2GHzB\ ' 
filename{20}  =  'Air  test  10MHz  2GHzB  02  25  14  06  35  53'; 

dir{21}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz__2GHzC\ ' 
filename{21}  =  'Air  test  10MHz  2GHzC  02  27  14  06  43  55'; 

dir{22}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz_2GHzd\ ' 
filename{22}  =  'Air  test  10MHz  2GHzd  03  06  14  08  52  04'; 

dir{23}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  3GHz\'; 
filename{23}  =  'Air  test  10MHz  3GHz  03  11  14  06  32  51'; 

dir{24}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz_3GHzB\ ' 
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filename { 24 } 


=  ' Air_test_10MHz_3GHzB_03_13_14_06_32_45 ' ; 
dir{25}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  3GHzC\ ' ; 
filename{25}  =  'Air  test  10MHz  3GHzC  03  18  14  06  32  08'; 

dir{26}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  3GHzd\ ' ; 
filename{26}  =  'Air  test  10MHz  3GHzd  03  20  14  08  09  07'; 

dir{27}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  4GHza\ ' ; 
filename{27}  =  'Air  test  10MHz  4GHza  04  10  14  07  01  22'; 

dir{28}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  4GHzb\ ' ; 
filename{28}  =  'Air  test  10MHz  4GHzb  04  11  14  06  10  18'; 

dir{29}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  4GHzc\ ' ; 
filename{29}  =  'Air  test  10MHz  4GHzc  04  15  14  07  01  09'; 

dir{30}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  10MHz  4GHzd\ ' ; 
filename{30}  =  'Air  test  10MHz  4GHzd  04  16  14  06  42  27'; 

dir{31}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air  test  100MHz  lGHz\ ' ; 
filename{31}  =  'Air  test  100MHz  1GHz  02  19  14  07  20  46'; 

str_meta  = 

sprintf ( ' %s%s .mspecrawdata ' , dir{SELECT} , filename { SELECT } ) 
str_data  = 

sprintf ( ' %s%s . specrawdata ' , dir { SELECT } , filename { SELECT } ) ; 

%  Get  metadata  on  selected  data  file 

fid  meta  =  fopen (str^meta) ; 

META  =  textscan(fid  meta, '%s' ) ; 

ave  =  META{ 1 } {1}; 

ref  =  META{ 1 } {2}; 

attn  =  META{ 1 } {3}; 

rbw  =  META{ 1 } {4}; 

%  Read  in  data  file 

fid_data  =  fopen (str_data) ; 
g=0; 
f =0  ; 
a=0  ; 

while (g==0 ) 

ft=fgetl (fid_data) ; 
f= [f ,  str2num (ft) ] ; 
at=fgetl (fid_data) ; 
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a= [a, str2num (at) ] ; 
g=feof (fid_data) ; 

end 

f=f (2 : end) ; 
a=a (2 : end) ; 
alength  =  length (a); 
display  (alength) ; 

%  Determine  resolution  bandwidth 

Bandwidth  =  f(3)  -  f(2); 

%  Plot  data  file 

figure  ( 1 ) ; 
plot (f /le6, a) ; 
axis  tight; 
grid; 

xlabel ( ' Frequency  (MHz )  ' )  ; 
ylabel ( ' Power  (dBm) ' ) ; 

title (strcat (rbw, ' , ave, ' , ref, 1 -- 
' , attn) ,  Interpreter  ,  none ' ) ; 


%  Modify  RF  spectrum  based  on  the  Watershed  Transform  concept 
%  Locate  peak  values  in  the  RF  spectrum 
%  Fill  in  valleys  between  the  peak  locations 


Inputlndex  =  alength; 
InputArray  =  a; 

Input  =  a; 

WatershedArray  =  InputArray; 


i  =  1; 

tempi  =  InputArray ( 1 ) ; 
negativelocation  =  0; 
peaklocationl  =  0; 
peaklocation2  =  0; 

while  i  <  Inputlndex-l 

while  i  <  Inputlndex-l 

if  InputArray ( i+1 )  ==  InputArray ( i ) 
i  =  i + 1  ; 

elseif  InputArray ( i+1 )  >  InputArray ( i ) 
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tempi  =  InputArray ( i+1 ) ; 
peaklocationl  =  i+1; 
i  =  i  +  1 ; 


else 

temp3  =  InputArray ( i+1 ) ; 
negativelocation  =  i+1; 
break 

end 


end 

while  i  <  Inputlndex-l 

if  InputArray ( i+1 )  <  InputArray ( i ) 
temp3  =  InputArray ( i+1 ) ; 
negativelocation  =  i+1; 
i  =  i  +  1  ; 


else 

k  =  i; 

m  =  negativelocation; 

i=  i  -  1; 

while  m  <  Inputlndex-l 

if  InputArray (m+1 )  ==  InputArray (m) 
m  =  m  +  1  ; 

elseif  InputArray (m+1 )  >  InputArray (m) 

temp2  =  InputArray (m+1 ) ; 
peaklocation2  =  m+1; 
m  =  m  +  1  ; 


else 

break 

end 


end 

break 

end 


end 

k  =  peaklocationl; 


if  abs (peaklocation2-peaklocationl )  ==  2 
if  tempi  <  temp2 

WatershedArray ( k+1 )  =  tempi; 


else 

WatershedArray ( k+1 )  =  temp2; 

end 
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else 


if  tempi  <  temp2 

while  InputArray ( k+1 )  <  tempi 

WatershedArray ( k+1 )  =  tempi; 
k  =  k+1; 

if  k  ==  Inputlndex 
break 

end 

end 

else 

while  InputArray ( k+1 )  >  temp2 

k  =  k+1; 

if  k  ==  Inputlndex 
break 

end 

end 


for  c  =  k+1 : peaklocation2 

WatershedArray (c)  =  temp2; 

end 

end 

end 

tempi  =  InputArray ( i+1 ) ; 
peaklocationl  =  i  +  1; 
i  =  i  +  1  ; 

end 

%  Determine  some  statistics  values  for  data 

a  =  WatershedArray; 

M  =  mean  (a) ; 

Med  =  median (a); 

S  =  std (a) ; 

Max  =  max (a) ; 

Min  =  min (a) ; 

Kurt  =  kurtosis (a) ; 

Range  =  abs (Min  -  Max) ; 

Number  =  floor (Range) ; 
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Bins  =  Number  *  2; 


disp([M  S  Med  Max  Min  Kurt] ) ; 


%  Initialize  parameters 

Erosionlndex  =  1; 

InputArray  =  a; 

ErosionArray  =  a; 

DilationArray  =  a; 

Medsize  =  5; 

Medcount  =  fix (Medsize) ; 

Medstart  =  Medcount  +  1; 

%  Process  data  file  with  a  running  median 

for  m  =  Medstart : Inputlndex-Medcount 

InputArray (m)  =  median (a (m-Medcount :m+Medcount) ) 

end 

for  m  =  l:Medcount 

InputArray (m)  =  InputArray (Medstart) ; 

end 

for  m  =  Inputlndex-Medcount : Inputlndex 

InputArray (m)  =  InputArray ( Inputlndex-Medcount)  ; 

end 


%  Plot  processed  median  of  data  file 

figure (2 )  ; 

plot (InputArray) ; 

PreviousPower  =  sum (a) /Inputlndex; 
display  (PreviousPower) ; 

%  Initialize  parameters 

j=  i; 

Diff Power  =  1; 

DilationMedian  =  Med; 

DilationMax  =  Max; 

DilationCF  =  DilationMax-DilationMedian; 
%  Perform  Erosion  operation  on  data  file 
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if  DilationCF  >  10  &&  j  ==  1 
display ('CF  &  Power'); 


while  DilationCF  >  10 

for  k  =  1 : Inputlndex- j 
ErosionArray ( k)  =  min ( InputArray ( k : k+ j ) ) ; 
end 


for  k  =  Inputlndex- j : Inputlndex 

ErosionArray ( k)  =  InputArray ( Inputlndex- j ) 

end 


%  Plot  results  of  Erosion  operation 

figure  ( 5 ) ; 

plot (ErosionArray) ; 


%  Perform  Dilation  operation  on  data  file 
DilationArray  =  ErosionArray; 
for  k  =  Inputlndex : j +1 

DilationArray ( k)  =  max (ErosionArray ( k-j : k) ) 

end 


for  k  =  j : 1 

DilationArray ( k)  =  ErosionArray ( j +1 ) ; 

end 

%  Plot  results  of  Dilation  operation 

InputArray  =  DilationArray; 
figure ( 6) ; 

plot  (DilationArray) ; 

DilationPower  =  sum (DilationArray) /Inputlndex; 
DilationMedian  =  median (DilationArray) ; 
DilationMax  =  max (DilationArray) ; 

DiffPower  =  abs ( PreviousPower  -  DilationPower); 
DilationCF  =  DilationMax-DilationMedian; 

Approved  for  public  release;  distribution  is  unlimited. 

23 


display  ( [ j  DilationPower  DilationMedian  DilationMax 
DilationCF  DiffPower] ) ; 

PreviousPower  =  DilationPower; 

j=  j+i; 

end 

%  Check  performance  of  the  Erosion  &  Dilation  to  see  if 
operations 

%  resulted  have  converged  to  accept  criteria  (If  power 
difference  in 

%  the  RF  spectrum  is  less  than  0.1  between  iterations  previous 
and 

%  current  processing  of  RF  spectrum 
while  DiffPower  >  0.01 


display ( ' Power ' ) ; 


for  k  =  1 : Inputlndex- j 

ErosionArray ( k)  =  min ( InputArray ( k : k+j )) ; 


end 

for  k  =  Inputlndex- j : Inputlndex 

ErosionArray ( k)  =  min ( InputArray ( k : Inputlndex) )  ; 

end 


figure  ( 5 ) ; 

plot (ErosionArray) ; 


DilationArray  =  ErosionArray; 

for  k  =  Inputlndex : j +1 

DilationArray ( k)  =  max (ErosionArray ( k-j ; k) ) ; 

end 

for  k  =  j : 1 

DilationArray ( k)  =  max (ErosionArray ( k : 1 )) ; 
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end 


InputArray  =  DilationArray; 
figure ( 6)  ; 

plot  (DilationArray) ; 

DilationPower  =  sum (DilationArray) /Inputlndex; 
DilationMedian  =  median (DilationArray) ; 

DilationMax  =  max (DilationArray)  ; 

DiffPower  =  abs (PreviousPower  -  DilationPower); 
DilationCF  =  DilationMax-DilationMedian; 

display  ( [ j  DilationPower  DilationMedian  DilationMax 
DilationCF  DiffPower] ) ; 

PreviousPower  =  DilationPower; 

j=  j+i; 

end 

else 


while  DiffPower  >  0.01 
display (' Power  Only'); 


for  k  =  1 :  Inputlndex- j 

ErosionArray ( k)  =  min ( InputArray ( k : k+j )) ; 

end 

for  k  =  Inputlndex- j : Inputlndex 

ErosionArray ( k)  =  min ( InputArray ( k : Inputlndex) ) 

end 


figure  ( 5 ) ; 

plot (ErosionArray) ; 


DilationArray  =  ErosionArray; 

for  k  =  Inputlndex : j +1 

DilationArray ( k)  =  max (ErosionArray ( k-j ; k) ) ; 

end 
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for  k  =  j : 1 

DilationArray ( k)  =  max (ErosionArray ( k : 1 ) ) ; 

end 

InputArray  =  DilationArray; 
figure ( 6 ) ; 

plot (DilationArray) ; 

DilationPower  =  sum (DilationArray) /Inputlndex; 
DilationMedian  =  median (DilationArray) ; 

DilationMax  =  max (DilationArray)  ; 

DiffPower  =  abs (PreviousPower  -  DilationPower); 
DilationCF  =  DilationMax-DilationMedian; 

display  ( [ j  DilationPower  DilationMedian  DilationMax 
DilationCF  DiffPower] ) ; 

PreviousPower  =  DilationPower; 

j=  j+i; 

end 


end 


%  Plot  original  RF  spectrum  and  threshold 

figure  ( 3 ) ; 

plot (f/le6, Input) ; 
axis  tight; 

grid; 

xlabel ( ' Frequency  (MHz )  '  )  ; 
ylabel ( ' Power  (dBm) ' ) ; 

title (strcat (filename (SELECT) , '  --  ' , rbw, ' 

Hz'),'  Interpreter ' ,  ' none ' ) ; 
hold  on 

%  Determine  threshold  array 

Threshold  =  DilationArray+25- (2 . 9*logl0 (Bandwidth)); 


plot  (f/le 6, Threshold,  ' r ' ) ; 


hold  off; 

%  Save  results  to  data  file 

lab  =  num2str (SELECT)  ; 
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label  =  strcat (' WatershedDetection-Figure ', lab, '. jpg ') ; 
saveas (gcf , label)  ; 

%  Plot  original  RF  spectrum  and  final  results  of  Erosion  & 
Dilation 
%  processing 

figure (10)  ; 

plot (f/le6, a) ; 
axis  tight; 

grid; 

xlabel ( ' Frequency  (MHz )  '  )  ; 
ylabel ( ' Power  (dBm) ' ) ; 

title (strcat (filename (SELECT) , '  --  ' , rbw, ' 

Hz'),'  Interpreter ' ,  ' none ' ) ; 
hold  on 


plot (f/le6, DilationArray,  ' — r ' )  ; 


hold  off; 


label  =  strcat (' WatershedDilation-Figure ', lab, '. jpg ') ; 
saveas (gcf, label) ; 

%  Plot  orginal  RF  spectrum  and  final  results  of  Erosion  and 
Dilation 
%  processing 
figure  (40) ; 


plot (f/le6, Input) ; 
axis  tight; 
grid; 

xlabel ( ' Frequency  (MHz )  ' )  ; 
ylabel ( ' Power  (dBm) ' ) ; 

title (strcat (rbw,  -  - ' , ave,  ' , ref,  -- 
' , attn) , ' Interpreter  ,  none ' ) ; 

hold  on 

plot (f /le6, WatershedArray, ' r ' ) ; 
hold  off; 


end 
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Intentionally  lelt  blank. 
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Appendix  B.  Graphs  of  Morphological  Processed  RF  Spectrum 

Files 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


ARL 

US  Army  Research  Laboratory 

CF 

crest  factor 

GHz 

gigahertz 

kHz 

kilohertz 

MHz 

megahertz 

PDF 

probability  density  function 

RBW 

resolution  bandwidth 

RF 

radio  frequency 
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